Simulation of Potts models with real q and no critical slowing down 
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A Monte Carlo algorithm is proposed to simulate ferromagnetic g-state Potts model for any real 
Q > 0. A single update is a random sequence of disordering and deterministic moves, one for each 
link of the lattice. A disordering move attributes a random value to the link, regardless of the state 
of the system, while in a deterministic move this value is a state function. The relative frequency 
of these moves depends on the two parameters q and (5 = The algorithm is not affected by 
critical slowing down and the dynamical critical exponent z is exactly vanishing. We simulate in 
this way a 3D Potts model in the range 2 < g < 3 for estimating the critical value qc where the 
thermal transition changes from second-order to first-order, and find qc = 2.620 ± 0.005 . 

PACS numbers: 05.50.-Fq, 64.60.Fr,75.10.Hk 
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I. INTRODUCTION 

Q-state Potts model is perhaps one of the simplest, 
non-trivial models in statistical mechanics. A broad set 
of techniques has been brought to bear on it in a variety 
of disciplines and it has been the subject of considerable 
theoretical attention over the last two decades (for a re- 
view, see dl). 

This model is theoretically well defined for any real or 
complex value of g In particular, the limit g — > 1"*" 
corresponds to the random percolation problem and the 
limit <7 — > 0"*" has a fundamental role in enumerating the 
spanning trees of a graph ||^. Two-dimensional confor- 
mal field theory Q suggests exact formulae for the critical 
indices and for other universal quantities as continuous 
functions of q in the range < g < 4. Another interest- 
ing problem involving non-integer q in three-dimensional 
Potts models is the determination of the universal value 
qc for which the thermal transition changes from second- 
order to first-order. A variety of techniques have been 
used 1^-^, which locate qc in the range 2 < < 3. All 
these methods require extrapolations in q because the 
standard simulations work only at integer values of q. 
Reweighting techniques ||ic| , [ll|] and transfer matrix meth- 
ods |12[| allow to estimate some thermodynamic functions 
1^ in a wider range of q, however there is no way to eval- 
uate correlation functions there. 

In this paper we remove this limitation by construct- 
ing a new Monte Carlo (MC) algorithm which works for 
any real g > 0. Although the time required for a sweep 
through the system grows faster than its size because 
at some step of the algorithm nonlocal information is 
required, the simulations are not affected by a critical 
slowing down and the dynamical critical exponent z is 
exactly zero. We test the reliability of the method by 
comparison with some exact results for 2D Potts model 
at criticality. We probe its effectiveness by performing 
large scale MC simulations of a three-dimensional Potts 
model for estimating the universal value qc- 



= 1,2,. 



II. THE ALGORITHM 

Starting with the Hamiltonian H = 
where the site variable Ui takes the values ui 
with [ij) ranging over the links of an arbitrary lattice or 
graph A, one can write the g-state Potts model partition 
function Z = X]{cr} in the Fortuin Kasteleyn (FK) 

random cluster representation H 



Y,W{G) = Y,n{b,c)v\^ 



(1) 
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where v = — 1 ^ summation is over all span- 

ning subgraphs G C A, W{G) — v^q'^ is their weight, 
expressed in terms of the number 6 of edges of G, called 
bonds, and the number c of connected components or FK 
clusters; 17(6, c) is the number of subgraphs with b bonds 
and c clusters. This representation now defines a model 
for any real or complex q. 

In principle, one could directly use Eq.(^) to define a 
Metropolis algorithm working for positive non integer q 
list , but this is a difficult problem to simulate because, 
for each proposed change of a link, the number c of FK 
clusters, a nonlocal property, must be determined. Large 
lattices require a huge amount of CPU time. As a matter 
of fact, such a method has been applied only to two- 
dimensional systems, where special topological relations 
can be used [ p^ . 

Our strategy is different. We start by considering a 
useful identity which can be derived using the methods 
described in Ref. Q. 

Let I be any link of A. Denote by {G^} the set of 
spanning subgraphs where / is a bond and by {G~[} those 
in which this bond is missing. We have Z = Zf -I- Z'[ , 
with Zf' — X^G* W^l^*;^)- Introducing a bond varialale 
ai equal to 1 when / is a bond and otherwise yields 



Z^ 



(2) 
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The same quantity can be evaluated in a different way 
by addition of a bond to each graph of type ■ There 
are two kinds of missing bonds. Those joining two differ- 
ent clusters, called potential bridges, are picked out by a 
variable f3i which takes the value 1 only on them and is 
zero otherwise; their addition lowers the number c of FK 
clusters. We have 

A = l^ W{G^)^^W{G+) . (3) 

The remaining missing bonds, described by a similar vari- 
able 7/, join two sites of the same cluster; their addition 
keeps c invariant, then 

7i = 1 ^ WiGi)v = WiG+) . (4) 

Combining Eq.s (|), (|) and (§ yields (a,) = + 
v{'yi) which is the wanted identity. Since of course ai + 
/?; + 7i = 1 7 we can rewrite it as 

{ai)=p{ai) + ^{Pi)+p{ji) , (5) 
Q 

where the weighting factors can now be interpreted in 
terms of probabilities. The idea is now to regard this 
identity as the limit of a recursion relation of the type 

7r("+i) =pa(") + ^/3(«)+p^(") , (6) 

where 7r["~''^'' is the probability of having a bond on the 
link I in the configuration G*^"+^^ It is expressed as 
a state function {ai,Pi, or 7/) of the same link in the 
G^"^ configuration. This generates a Markov process 
•• — > G*'"-' G*^""'"^-' ■■ where the equilibrium distri- 
bution yields Eq. (^. This stochastic chain fullfills two 
important conditions: i) there is a non-zero probability 
of going from any configuration to any configuration in a 
single sweep through A, ii) the equilibrium distribution 
maps to itself as Eq.(||) is kept invariant by the process. 
One can then argue that detailed balance is satisfied. 

To see it directly, assume for instance that in the n*'* 
configuration / is a potential bridge = 1 , G^"^ = 

Gf) which is promoted to a bond in the (n -I- 1)*'* con- 
figuration (a["+^^ = 1 , G("+i) = G+). The transition 
rate is P{G^ Gf) = ^. Conversely Eq.(|) yields 
P{Gj G;") = 1 - p. Then, according to Eq.s (0) and 

P{Gi^Gt) ^ W{Gt) 
P{Gt^Gi) W{Gi) 

as detailed balance requires. The same conclusion can be 
reached in all the other cases. 

A straightforward, preliminary, implementation of the 
recursion relation (||) is the following: i) go over each link 
Z G A of the configuration G^"-* and generate a pseudo- 
random number Xi uniformly distributed from to 1. ii) 



Create a bond on / only in the following two cases: a) 
Xi < p and I is a bond {ai = 1) or a missing bond joining 
two sites of the same FK cluster (7; = 1) ; b) X; < | and 
I is a potential bridge (/?/ — 1). This generates uniquely 
the configuration g'"^^'. 

Let g > 1 for definiteness. It is worth noting that when 
AT/ < I the algorithm adds a bond to I regardless of 

which configuration G*-"' we are dealing with. Similarly, 
when Xi > p no bond is added. In the remaining cases 
< Xi < p). The value attributed to I (bond or no 

bond) is unambiguously determined by G^"'' . 

Inspecting all the cases leads to the following better 
implementation of the algorithm: 

Step 1: Pick a link / e A and generate a pseudo- 
random number Q < Xi < 1. 

Step 2: Update the link according to the following 
scheme 

move current state new state 
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Step 3: Return to step 1. 

The first two moves do not need any information on 
the state of the system: they just disorder it. The last 
one is a purely deterministic move; its only effect is to 
put a bond whenever a link joins two sites of the same 
cluster. It requires distinguishing between the two kinds 
of missing bonds {Pi = 1 or 7; = 1). One can infer this 
nonlocal property by identifying the connected compo- 
nents of the configuration, like in the Swendsen-Wang 
(SW) algorithm |jl^. This cluster reconstruction is time 
demanding, however it gives a complete information on 
the state of the missing bonds of the whole lattice. As 
the update proceeds through the lattice this amount of 
information is progressively lost because of disordering 
moves (the deterministic moves never change c) . We may 
partly keep track of the cluster structure by relabelling 
the cluster indices whenever a disordering move creates 
a bond between two of them. Cluster reconstruction is 
truly necessary only when a deterministic move touches 
a missing bond of a putative single cluster where some 
bond has been erased by previous disordering moves. 

Because of non locality, the number of operations in- 
volved every MC step is oc N", where TV is the number 
of links and 1 < a < 2. The efficiency of the algorithm 
depends crucially on the actual number of cluster recon- 
structions per sweep. In our 3D simulations reported be- 
low the fraction of links requiring cluster reconstruction 
was about 3% with a decreasing trend for larger lattices. 
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TABLE I. The decorrelation time r of the new algorithm 
for the critical 2D Ising model for different linear lattice sizes 
L is compared with the same quantity of the SW algorithm 
Tsw and with the upper bound To- The definition of r and 
the TSW data are taken from Ref |16|. 



L 


T 


Tsw 


To 


8 


2.65(3) 


5.17696(32) 


3.3869 


16 


3.16(5) 


6.5165(12) 


4.5158 


32 


3.69(6) 


8.0610(18) 


5.6448 


64 


4.3(1) 


9.794(4) 


6.7737 
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III. CORRELATION TIMES 

An unusual feature of the described algorithm is the 
presence of randomly distributed disordering moves. The 
mean number of links subjected to disordering moves in 
a single sweep is Npr with pr = 1 + ^ — p . For in- 
stance, in the Ising model {q = 2) at criticality more 
than 70% of the links are disordered every sweep. It is 
now easy to find an upper bound for the mean number r 
of MC steps needed to generate effectively independent 
configurations. After n sweeps the mean number of links 
which do not have yet undergone a disordering move is 
N(l — Pr)". When this number is of the order of 1 all 
the links have been touched by a disordering move and 
the upper bound To > t is given by the obvious relation 
Nil-Pry" ~ 1, i.e. To = - j,g'^f5p„) - Thus the dynam- 
ical exponent 2; is 0, as critical slowing down manifests 
itself by the power law r cx N'^ at the critical tempera- 
ture where a second order phase transition occurs 
A numerical estimate of the decorrelation time of the dy- 
namics of the new algorithm for the critical Ising model 
on a square lattice is reported in Tab. |. Note that the 
actual value of r does not saturate the upper bound and 
is much smaller than the analogous quantity of the SW 
algorithm. 

The new algorithm proves also useful in fighting 
against another dynamical problem which one deals with 
in the case of first-order transitions, namely the ex- 
ponentially fast suppression of the tunnelling between 
metastable states with increasing lattice size. 

To reduce this type of slowing down the multi- 
canonical MC algorithm has been proposed ; also the 
method of simulated tempering proves useful |19|. In 
a few numerical tests for two-dimensional models with 
q — 7 and 20 > L > 100 we found that the tunnelling 
time of the canonical algorithm described in the present 
paper grows with the system size as oc with 
a = 1.03 ± 0.03, like in an optimal variant |^ of the 
multi-canonical method, but with a smaller proportion- 
ality factor. 

The reason of this performance is that the random 
moves accelerate the tunnelling between order and disor- 
der. The drawback is that the new algorithm is non-local. 




FIG. 1. Plot of A{E,L) resulting from a simulation of a 
3D Potts model at g = 2.75 with 3.3 lO'' MCS for L = 14 (full 
circles) compared with the extrapolation at the same value 
of q of an actual SW simulation of 4.1 lO'' MCS at g = 3 
(crosses). The latter data are shifted to the right for clarity. 

so the CPU time grows like with 6 > 1; for instance 
in the present case we found b ^ 1.85. Thus the new 
algorithm cannot certainly be recommended for integer 
q, although at a first-order transition it performs much 
better than any local canonical algorithm. 



IV. SIMULATIONS 

As a first, simple, application of the new algorithm we 
tested the reliability of our code by checking a percola- 
tion property of Potts model on a square lattice which 
is supposed to be exact in the range < g < 4, namely 
that the mean frequency of active bonds (ai) at critical- 
ity (corresponding to v = y/q in Eq.(|^)) should coincide, 
in the thermodynamic limit, with the random percola- 
tion value, i. e. (a;) — ^, irrespective of the value of q 
We simulated critical Potts models in a 128 x 128 
square lattice with q ranging from 1.5 to 3.5. In all the 
cases the mean number of bonds was compatible with the 
exact result. Finite size effects of this observable, which 
are visible on smaller lattices, allow to evaluate the criti- 
cal thermal exponent v as a function of q. This could be 
used to check a conjectural formula suggested by the 2D 
conformal field theory Q . We plan to study this problem 
in a future publication. 

The irew algorithm allows us to deal with an important 
issue of three-dimensional Potts model, namely the esti- 
mate of the tricritical point qc in the range 2 < < 3, 
where the thermal transition changes from second-order 
for q < qc to first-order for q > q^. Many different tech- 
niques have been used to locate this point (5|-||]. We ap- 
plied a method very similar to that described by Lee and 
Kosterlitz jsj by computing the double histogram N{h, c) 
of bond and cluster number distribution in a cubic lattice 
of volume at a given (3 and q and then extrapolating 
the data to nearby values. Using Eq.(|^) we can write 
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TABLE IL The simulations were performed on cubic lat- 
tices of side L at the values of q and /3 listed below. MCS is 
the number of Monte Carlo steps considered. 



0.004 



L 


1 


13 


MCS 


0.003 


12 


2.70 


0.52270 


3.0 10' 


P{E) 


13 


2.70 


0.52270 


3.0 lO'' 


14 


2.75 


0.52721 


3.3 lO'^ 


0.002 



N{h,c-(3,q,L)=Un{b,c) 



(8) 



where N is the number of MC sweeps. We can trade 
the number of bonds b for the energy per site E using 
the relation E ~ — fo^^. Near a first-order transition 
the histogram P{E) — N{h, c)/J\f has a characteristic 
double peak structure corresponding to the ordered and 
the disordered phase. A suitable reweighting through 
Eq.(||) of the energy distribution yields the value (3c{L, q) 
where the two peaks at Ei {(3, L) and E2 {(3, L) are of equal 
height. A typical plot of the quantity A{E,q;l3c,L) = 
— ln(iV(6, c) /Af) is shown in Fig.l. A useful estimator 
of the interface free energy between the ordered and the 
disordered phase [|T] is given by 



AF{q, L) = A{E,r^,q■, Pc, L) - A{Euq; (3c, L) 



(9) 



where Em is the local maximum which separates the two 
dips at El and E2 (see Fig.l). At a first-order transition 
AF{L) increases monotonically with L and is expected 
to vanish at the tricritical point. Extrapolating the nu- 
merical data both in (3 and q one may locate this point. 
The region of reliable extrapolation |11| is 0{1/L^) for 
both /3 and q. This does not cause a problem for since 
it can be adjusted continuously, but q cannot in standard 
simulations, being by necessity an integer value. Actu- 
ally Lee and Kosterlitz performed their simulations at 
q = 3 and found that the extrapolated data become too 
noisy for \6q\ > 0.3 ||2^. In our case we can directly 
evaluate the range of reliable extrapolations. Indeed the 
main advantage of the algorithm described in this paper 
is that now also q can be adjusted continuously. 

Our simulations were performed on three different lat- 
tices as listed in Tab.^. The statistics is good since in 
all the cases the mean flipping time between coexisting 
states was no larger than thirty MC steps. The errors 
were calculated by gathering the histogram N(b, c) every 
10^ MC steps and then performing a standard analysis. 

In all the cases the energy histogram showed a double 
peak structure, providing a direct evidence of the first- 
order nature of the transition for these values of q (see 
Fig. 2). This yields the upper bound qc < 2.7. In shorter 
simulations at q = 2.6 we found no trace of a double peak 
structure. This suggests qc > 2.6. Using the reweighting 
method we estimated the values (3c{L,q) where the two 
peaks are of equal height for each L and for few values of 




0.001 



FIG. 2. Energy histogram at f3c{L,q) obtained by a sim- 
ulation at g = 2.7 and L = 13. 

q near q = 2.7 and the corresponding values of AF. The 
results are reported in Fig. 3. 

A further reweighting up to (7 = 3 allowed to compare 
the extrapolated data with those coming from a similar 
extrapolation of standard SW simulations ed q — 3. This 
comparison showed that the range of reliable extrapo- 
lations is \Sq\ < 0.25. It has to be noted that we could 
not use for this comparison the high precision data of Ref. 
[ p5[ , because the energy distribution in terms of spin vari- 
ables used there does not coincide with that expressed in 
terms of bond variables used by necessity in the present 
approach. In particular the PdL, (?)'s are shifted and our 
AF{q, L) is always smaller. 

Simple finite size scaling considerations suggest § 
that near qc the interface free energy has the simple 
form AF{q,L) ~ (9 — qc)^L°' which fits very well to 
our data (see Fig. 3). To within our numerical accuracy 
a = 4.8±0.1 and qc = 2.620±0.005. This agrees with the 
value qc = 2.55±0.12 obtained in the large q expansion of 
the latent heat |6). Lee and Kosterlitz |^ extrapolating 
q = 3 data found a smaller value, qc — 2.45 ± .0.10. The 
difference could be due to the fact that extrapolations 
with \5q\ > 0.25 gives an overestimate of AF (this is al- 
ready visible in Fig.l). Other approximate methods give 
even smaller values: real space renormalization group 
methods [|j yield qc ~ 2.2 while an Ornstein-Zcrnike ap- 
proximation M gives qc ^ 2.15. 



V. CONCLUSIONS 

This work provides a new MC algorithm to simulate 
ferromagnetic q— state Potts model which has two very 
unusual features: it works for any real q > and does 
not suffer of any critical slowing down. The former prop- 
erty is an obvious consequence of the fact that it is based 
on the Fortuin Kasteleyn random cluster representation, 
where q acts as a continuous parameter. The latter is 
more tricky and is due to the implementation of the al- 
gorithm with a random sequence of disordering moves. 
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2.7 2.7.5 2.8 2.85 

<1 

FIG. 3. Plot of AF((7, L) near q = 2.7 . 

randomly distributed over the lattice. There is no rea- 
son to believe that this disordering mechanism is specific 
to Potts model and it would be very interesting to try 
to implement it in other, more general MC methods. A 
drawback of the new algorithm is that it is non-local, so 
the CPU time of a single sweep grows with the volume 
V as with 1 < 6 < 2, thus it is not recommended 
for integer q, where the SW algorithm works with 6=1. 
Actually at a first-order transition the new algorithm per- 
forms better than the SW method, but there the multi 
canonical MC algorithms are more suitable. 

It is straightforward to extend the new algorithm in 
order to take into account quenched bond randomness, 
provided that all the couplings are ferromagnetic. On the 
contrary, generalising to systems with frustrations seems 
a rather difficult task, because it is not obvious how to 
define in this case the FK clusters for non- integer q . 

We used such an algorithm to study the region 2 < 
g < 3 of a three-dimensional Potts model in order to esti- 
mate the critical value qc for which the thermal transition 
changes from second to first-order. We obtain a rather 
precise estimate compared to other methods, the 
reason being that all the other methods are based on 
extrapolations from integer values of q, while the new 
algorithm simulates the system at nearby values of qc- 

The author would like to thank M. Caselle and A. 
Coniglio for helpful discussions. 
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